Overlapping resonances in the control of intramolecular vibrational redistribution 
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Coherent control of bound state processes via the interfering overlapping resonances scenario 
[Christopher et al, J. Chem. Phys. 123, 064313 (2006)] is developed to control intramolecular 
vibrational redistribution (IVR). The approach is applied to the flow of population between bonds 
in a model of chaotic OCS vibrational dynamics, showing the ability to significantly alter the extent 
and rate of IVR by varying quantum interference contributions. 

PACS numbers: 



INTRODUCTION 



Quantum control of molecular processes^ 2 - has proved, 
over the past two decades, to be viable both theoretically 
and experimentally. An examination of the coherent con- 
trol literature, wherein scenarios are expressly designed 
to take advantage of quantum interference phenomena, 
shows that the vast majority of applications has been to 
processes occurring in the continuum energy regime. Re- 
cently we proposed a new approach to controlling bound 
state dynamics in large polyatomic molecules^ that ex- 
ploits interferences between overlapping resonances. We 
have demonstrated the viability of this scenario in con- 
trolling internal conversion in pyrazine^^i In the present 
paper we further develop this method, applying it to 
the control of Intramolecular Vibrational Redistribution 
(IVR). As an example, we study the control of the flow of 
energy between bonds in a model of OCS. This molecule, 
though small, is of particular interest at high energies, 
where, classically, it displays predominantly chaotic dy- 
namics. In spite of the classical chaos, quantum control 
via the present scenario is shown to be excellent. 

This paper is organized as follows: Section II provides 
an overview of the theory, with a discussion of the Fes- 
hbach partitioning technique which, as we have shown, 3 
provides a highly efficient method for dealing with bound 
state problems. Section III describes the collinear OCS 
model and its classical dynamical characteristics. In Sec- 
tion IV we discuss the application of the method to the 
control of IVR in OCS. An Appendix describes our use of 
the Feshbach partitioning technique for the numerical so- 
lution of the bound state problem for small systems such 
as OCS. A more ambitious method for addressing consid- 



* Current address: Institute- de Matematicas y Fi'sica Fundamental, 
Consejo Superior de Investigaciones Cienti'ficas, Serrano 123, 28006 
Madrid, Spain 



erably larger systems, the "QP algorithm", is described 
elsewhere £ 



II. 



BOUND STATE CONTROL 



Time-evolution of populations in molecular 

systems 



We consider a system described by an Hamiltonian H 
which can be partitioned physically into the sum of two 
components Ha and Hb, plus the interaction Hab be- 
tween them: 



H = H A + H E 



H 



AB, 



(1) 



The eigenstates and eigenvalues of the full Hamiltonian 
are defined by: 



H\i) = £ 7 | 7 > 



(2) 



The ( "zeroth-order" ) eigenstates and eigenvalues of the 
sum of the decoupled Hamiltonians are defined as 



{H a + H b )\k) = E^\k) 



(3) 



Below, we are interested in the time evolution of the sys- 
tem, initially prepared in a superposition of zeroth order 
states. 



|MK* = o)) = 5> K | K >, 



(4) 



where {c K } are "preparation" coefficients. All sums over 
| At), here and below, are assumed to be confined to a 
subspace S. For example, the selected initial states might 
consist of a set with population heavily concentrated in 
one bond of a molecule, in which case, energy flow out of 
such superposition states is examined. 

The time-evolution of Eq. ((U at any subsequent time 
can then be obtained by expanding the (zeroth order) 



2 



eigenstates, |k), in terms of the exact eigenstates I7) to 
give: 

|*(t))=^ CK < 7e -^ t/,l |7), (5) 

with a* = (7|k). The structure of |(7|k)| 2 as a function 
of 7 defines a resonance shape that provides insight, in 
the frequency domain, into the population flow out and 
into the zeroth order \k) states. 

Given this time evolution, the amplitude for finding 
the system in a state \k) at time t is 

c K = («|*(t)> = CK'M K , K >(t), (6) 

k' 

where 

M K , K ,(t) = XX^e-^'/ft 

7 

= <*l ^E e ^* A l7)(7lj l«'> (7) 

is the (k, k') element of the overlap matrix M(t) defined 
by the term in brackets in Eq. |J7J|. Note that, for k' ^ K, 
if the states \k) and \k') do not overlap with a com- 
mon I7), i.e., there are no overlapping resonances, then 
M K>K i = 0. Our previous studies^ have demonstrated the 
significance of such overlapping resonances to the control 
of radiationless transitions, such as internal conversion. 

From Eq. ([6]), the probability of finding the system in 
a collection of states \k) contained in the initial set S at 
time t is given by 

F(t) = Xl<^l*W)| 2 =c t K(0c, (8) 

K 

where c is a K-dimensional vector whose components are 
the c K coefficients, and K(t) = M t (t)M(t). The gen- 
eralization to the question of finding population in an 
alternative collection of states, other than S, is straight- 
forward. However, it is unnecessary for the study below, 
as will become evident. Equation ([5]) allows us to ad- 
dress the question of enhancing or restricting the flow of 
probability out of S by finding the optimal combination 
of c K that achieves this goal at a specified time T. Ex- 
perimentally, the resultant required superposition state 
can be prepared using modern pulse shaping techniques. 



B. The Feshbach partitioning technique 

Our interest is to control the flow of population out of 
some generic molecular subspace into the entire molec- 
ular Hubert space. In order to do so we make use 
of the bound state version of the Feshbach partition- 
ing techniquei&i Here, since the control approach is be- 
ing tested on a small system, we solve the resulting 



equations in a straightforward way, as described in Ap- 
pendix A. Larger systems can take advantage of the "QP 
algorithm" .— 

The Feshbach partitioning technique is based on defin- 
ing two projection operators 

K 

which satisfy the following properties: 

Q 2 = Q, P 2 = P, (10a) 

[Q,P] = o, (10b) 

P + Q = I, (10c) 

where I is the identity operator. In what follows, the 
flow of probability of interest is from the Q space to the 
P space. 

Using Eqs. (|10c|l and ([9]), the eigenstates of the full 
Hamiltonian can be written as 

Similarly, the Schrodingcr equation can be expressed as 

[E 7 - H][P + Q]| 7 > = 0, (12) 

whereby multiplying it by P and then by Q, and using 
Eq. (fT0|) . one obtains the following set of coupled equa- 
tions: 

[£J 7 - PHP] P| 7 ) = PHQ\i), (13a) 
[E 1 -QHQ]Q\ 1 ) = QHP\i). (13b) 



The states \k) and |/3) are solutions to the decoupled 
(homogeneous) equations arising from Eqs. (|13ap and 
(|13b|) . respectively. That is, 

[Ep - PHP] P\P) = 0, (14a) 
[E k -QHQ]Q\k) =0. (14b) 

Contrary to continuum problems, in general £L =/= Ep 
and it is possible to express PI7) in terms of the partic- 
ular solution of the (inhomogeneous) Eq. (|13a|) . 

P| 7 ) = [P 7 - PHP]- l PHQYi). (15) 

Substituting Eq. (H} into Eq. (|13b|l results in 
[E 1 ~QHQ]Q\i}) =QHP[E^-PHP]- 1 PHQ\ip). (16) 
By rearranging terms in this equation, one obtains 

[E 1 -H]Q\ 1 ) = Q, (17) 
where TL is an effective Hamiltonian, defined as 

H = QHQ + QHP[E 7 -PHP]- 1 PHQ. (18) 
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i 


A 


Pi 


R° 


1 


0.08518 


1.5000 


2.9759 


2 


0.21238 


1.6251 


2.2559 


3 


0.16000 


1.1589 


2.8037 



TABLE I: Parameters defining the potential energy surface 
given by Eq. (|32[) . All magnitudes are given in a.u. 



The term between squared brackets can be written as 

[Ry - PHP}- 1 = J2 ^^HWI 
' ii 7 — tip 



by using the spectral resolution of an operator. The ma- 
trix elements of H. are given by 



(k|W|k'> = EJ K 



where 



-y 



E~ — En 



r K>lt > = 2irV(K\P)V(P\K'), 



(20) 

(21a) 
(21b) 



with V{k\(3) = (k\QHP\P) being the coupling term. 
Equations (|21a[) and (|21b[) represent the energy shift and 
the decay rate, respectively. By diagonalizing Eq. (|17p in 
a self-consistent manner, one obtains the energy eigen- 
values, E~f, and the values for the overlap integrals, a Kj7 . 

Note that the energy eigenvalues and the overlap in- 
tegrals can also be obtained^ by directly diagonalizing 
the full Schrodinger equation in the zeroth-order basis. 
However, the partitioning technique presented above has 
computational advantages for cases where the dimensions 
of the P space is large, since one only needs to diagonal- 
ize an effective Hamiltonian TL with dimensions given by 
the Q space. Note, however, that diagonalizing H re- 
quires using iterative procedures, and needs to be solved 
repeatedly until each eigenvalue is found. Appendix |A] 
provides details on the partitioning algorithm used here. 



C. Optimal Control 



To determine the set of optimal preparation coefficients 
leading to either a maximum or a minimum population at 
time t — T, we need to find the cxtrema of the function 9 



P x {t) = c f K(t)c - Ac f c 



(22) 



with respect to the coefficients c, where A is a Lagrange 
multiplier added to assure normalization, i.e., 



E 



= i. 



(23) 



The optimum vector, ct, is obtained by differentiating 
Eq. (|22| with respect to the components of , and equat- 
ing the resulting expression to zero at time T, i.e., 



dPx(t) 



dc* 



t=T 



K«/,re(T)c« - Ac« 



0. 



(24) 



The optimum vector resulting from this procedure is a so- 
lution to the eigenvalue problem represented by Eq. (|24| . 
Note that this vector can either maximize or minimize 
the solution. In the first case, the interference between 
overlapping resonances created by the initial superposi- 
tion will be seen to result in a delay of the population 
decay, whereas in the second case the population decay 
is being accelerated. 

In order to further clarify the dependence of the time- 
evolution on overlapping resonances, Eq. ([8]) can be re- 
expressed as 



K K. , K 



where 



and 



g K =^2\M K , iK \ 



(25) 



(26) 



(27) 



As expressed in Eq. (|25|) , the Q space population assumes 
the generic coherent control formi^: it is given as the 
sum of non-interfering pathways, represented by g K >, and 
interfering pathways, represented by f K ^ K . 



D. The Role of Overlapping Resonances 



The interference term in Eq. J25|) depends on f K ' tR , 
which, in accord with Eq. (|27[) . depends upon the over- 
lap between resonances. Qualitatively speaking, a reso- 
nance is described by bound states \k) coupled to a quasi- 
continuum of exact eigenstates I7). Each such state is 
thus associated with the energy width of the (k|7) over- 
lap coefficients. Overlapping resonances are the result of 
having at least two states whose resonance widths are 
wider than their associated level spacing. The resulting 
resonances interfere with one another, displaying a vari- 
ety of lineshapes?2 and are responsible for the interference 
in this control scenario. In the absence of overlapping 
resonances the full / K ' jt; -term in Eq. (|25[) vanishes and 
control disappears. 

Note that there are also contributions from overlap- 
ping resonances to the (? K -term, as can be seen from their 
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effect on the nature of the decay from the individual \k). 
These resonances distort the lineshape, and hence the 
corresponding time dependence. In order to determine 
the contribution from overlapping resonances, we have 
devised^ a qualitative measure, defined as 



P(t) = [P(t)-W(t)}, 



where 



w(t) = j2\ 



,M K . 



(28) 



(29) 



Here W(t) is a measure of the direct contribution, and 
P(t) provides a measure of the overlapping resonance 
contribution. In the absence of overlapping resonances, 
P(t) = W(t). 



III. CLASSICAL ASPECTS OF THE 
COLLINEAR OCS 



A. The OCS model 



As a working model to illustrate the usefulness of the 
method described in Sec. [ill we consider a collinear model 
of OCS, with a modified Sorbie-Murrell 10 potential. The 
interest in this system arises from the fact that, close to 
dissociation, i.e. in the energy region of interest below, 
the classical dynamics becomes highly chaotic. As such, 
collinear OCS is a complex system with a penchant for 
extensive IVR. 

The classical dynamics of OCS has been studied in 
both planar^ and collinea r 11 ' 12 versions. Here, we con- 
sider the collinear case, which is described by the Hamil- 
tonian 



H = 



p 2 



where 



2/X13 2/X23 

A*13 
A*23 



PlP2 

mc 



V{R U R 2 ,R 3 ), (30) 



m mc 

m + m c 
m s mc 



(31) 



are reduced masses; R\ and R2 are the CS and CO bond 
distances, respectively (R3 = R\ + R2)', and Pi and Pi 
are the corresponding momenta. 

In the course of this work we found that the Sorbie- 
Murrell OCS model 1 - displayed a second minimum at 
large distances along both the CS and CO exit channels. 
Although the depth of this second well is extremely small, 
there are a large number of closely packed eigenstates lo- 
calized in this region due to the length of the well. To our 
knowledge, there is no experimental evidence to either 
support or refute a second minimum, although they have 




FIG. 1: Contour plot of the potential energy surface given 
by Eq. (132[) . Solid and dashed lines represent, respectively, 
energy contours above and below the dissociation onset, at 
E = 0.100 a.u. (thick solid line). 



been associated^ with van der Waals interactions in O3. 
However, in order to ensure that the observed control is 
not a manifestation of this secondary minimum (as was 
the case in our preliminary studies) a modified interac- 
tion potential is used that removes these second minima 
while retaining the general features of the remaining po- 
tential. Specifically, the potential used here consists of a 
sum of three Morse functions, 



VXRi,i? 2 ,i? 3 ) 



3 3 

i=l i=l 



(32) 

with parameters given in Table |TJ A contour plot of the 
resultant potential energy surface is shown in Fig.[TJ Ex- 
cept for the Morse function V3, which depends on i? 3 , 
the parameters defining the other two Morse functions 
have been changed so that the potential smoothly fits 
the original one along the equilibrium directions while, 
at the same time, eliminating the second potential min- 
ima. Moreover, we have also modified the added con- 
stant in the potential so that the CS dissociation onset 
[V(oo, i?§) = D\ + D3] corresponds to the original value 
of E d = 0.100 a.u. 



B. Characterizing the Dynamics 

The classical dynamics of the resultant OCS model 
is characterized by a smooth transition from regular to 
chaotic dynamics with increasing energy. At an energy 
just below dissociation, (Eq = 0.097964 a.u., of interest 
below), the Poincare surface of section^ shows [Fig.[^a)] 
highly chaotic dynamics, with a stable region constitut- 
ing about 1/3 of the phase-space portrait. This energy 
corresponds to the mean value of the energies of the two 
wave packets, l^i), obtained below in maximizing and 
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in phase space grow in time: 



0.5 1 .0 

t(ps) 



A, 



t^oo t Uq 



(33) 



Here, in order to show how different regular and chaotic 
trajectories behave, we have computed the quantity 



A(t) = In 



d(t) 



(34) 



with do — 10 8 a.u. We label the finite time Lyapunov 
exponent, computed in this way, as At. 

The quantity A(t) is shown in Fig. HJb) for two sets 
of nearby trajectories^ picked in two different regions 
of phase-space: the stable island, and the chaotic sea. 



The results, to t « 1.2 ps, give A 



(stable) 



A 



(chaotic) 



1.46 ps , and 



17.41 ps 1 in the regular and chaotic regions, 



respectively. The associated times are to be compared 
to zeroth order vibrational periods (27.45 fs for the CS 
bond, and 18.10 fs for the CO bond). 

Finally, in Fig.[2jc) we show the energy in the CS bond, 
for a trajectory in the chaotic sea. As can be seen, the 
energy displays a complicated pattern, with irregular en- 
ergy transfer between both bonds as a function of time. 
Nonetheless, when one computes the energy average of 
an ensemble of trajectories, the pattern becomes smooth 
and displaying a profile than can be fitted to an expo- 
nential decays similar to those observed in its quantum 
counterpart below. 



FIG. 2: (a) Poincare surface of section for the collinear OCS 
model at E = 0.09796 a.u. The solid line represents the total 
energy contour, (b) Distance between two nearby trajecto- 
ries (with do = 10 -8 a.u.) chosen in the stable region (T), 
and in the chaotic sea (a). The high-frequency oscillations 
have been averaged out in both cases (the smoothing causes 

A (T) 

to appear as if it does not begin at zero), (c) CS-bond 
vibrational energy corresponding to the chaotic trajectory of 
part (b). 



minimizing the energy flow from the CS bond. Surfaces 
of section in the nearby energies are essentially similar. 
This being the case, there is no obvious classical origin 
to the control of bond energy relaxation described below. 
Of some future interest, however, might be an auxiliary 
study of the relationship of overlapping resonances in- 
duced control, observed below, to classical features such 
as bond energy recurrences, cantori, and the inhomoge- 
neous character of the OCS phase spac o 11 ! 12 ' 15 ! 16 ' 17 ! 18 . 

Quantitative insight into the rate of loss of correlations 
in the chaotic region of phase space can be obtained by 
computing Lyapunov exponents^ approximated by the 
average (over various trajectories) of the exponential rate 
at which the distance d(t) between adjacent trajectories 



IV. COHERENT CONTROL OF IVR 

A. Population decay 

We now consider the suppression (and enhancement) 
of IVR in the above model of OCS. Our intent is to assess 
the extent of control in such a system, and to establish 
the relationship between control and overlapping reso- 
nances. The coupling terms V(k\/3) and, subsequently, 
the overlap integrals a K . 7 and the energy eigenvalues E 7 
are calculated by expanding the OCS wave functions in 
products of the zeroth order states, 

m=T,^cs)®\^co)d mn ■ (35) 

m,n 

where \rfcs) anc ^ l£co) are eigenstates of the uncoupled 
CS and CO bond potentials, respectively, with quantum 
numbers m and n. Our interest is in the flow, for ex- 
ample, out of the CS bond. Hence, the Q subspace is 
chosen to represent all wave functions containing only 
excitation in the CS bond, i.e., \k) are \r]™ s ) ® |£co)' 
for all m, whereas the P subspace spans the space repre- 
sented by all other zeroth order excitations, i.e., the \j3) 
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IVR suppression IVR enhancement 
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0.01257 


- 0.08349 


0.00713 


20.76 


5 


0.0841783 


0.18017 


- 0.24251 


0.09127 


- 0.03560 


0.05674 


0.00449 


13.69 


6 


0.0837303 


0.20267 


0.15178 


0.06411 


0.19804 


0.05108 


0.04183 


34.02 


7 


0.0831998 


- 0.21171 


0.16534 


0.07216 


0.12120 


- 0.63185 


0.41392 


20.53 


8 


0.0825867 


0.25004 


- 0.41477 


0.23455 


0.20859 


- 0.39177 


0.19699 


25.32 


9 


0.0818910 


- 0.05482 


0.25131 


0.06616 


0.24895 


- 0.31440 


0.16082 


22.95 



TABLE II: Values corresponding to the eigenstates for the (uncoupled) CS bond used in the optimized superpositions. E K 
denotes the energy associated to these eigenstates; c r K and c\ are the real and imaginary parts of the c K coefficients, respectively; 
and tx is the decay time (see text for details). The optimization to maximize/minimize the energy transfer into the CO bond 
(suppression/enhancement of IVR) has been carried out at T = 100 fs. The energy corresponding to the ground state in the 
(uncoupled) CO bond is E° co = 0.00360475 a.u. 



are \Vcs) ® I Ceo)' n 7^ 0' describing excitation in the CS 
bond. Initiating excitation within Q and watching the 
flow into P then corresponds to an experiment wherein 
excitation flows out of the CS bond. 

As seen in Sec. MI A| the coupling term, QHP, nec- 
essary to obtain the energy shifts and decay rates, con- 
sists of a static term (V3), and a dynamic term [propor- 
tional to P1P2 in Eq. ([50)) ]. The overlap integrals and 
energy eigenvalues are obtained by self-consistent diago- 
nalization of Eq. (fl7|) . All vibrational states, \r)™ s ) and 
\£,qq), are numerically calculated using a discrete vari- 
able representation (DVR) technique, 21 obtaining a total 
of 45 eigenvectors for the CS bond, and 59 for the CO 
bond. The number of eigenvectors is larger in the second 
case, because the dissociation threshold of the CO bond 
is higher in energy. 

From all the vibrational states obtained, we have ob- 
served that control is best when considering a superpo- 
sition of states, i.e., Eq. (g]), that is near the dissocia- 
tion onset. The energy differences between these states 
are relatively small (« 0.0004 a.u., whose inverse cor- 
responds to a timescale of 60 fs), thus giving rise to 
a high density of states with time scales comparable to 
vibrational relaxation. The result is a greater opportu- 
nity for overlapping resonances which, as will be seen 
below, enhances the ability to control energy flow. In 
our case, the states used are the last nine bound eigen- 
vectors (before the dissociation onset) of the CS bond, 
whose corresponding eigenvalues are given in Table [TTJ 
Note, however, that dense eigenstate manifolds will oc- 
cur at far lower energies in larger molecules. Hence, the 
initial |$(0)) is comprised of a superposition of nine CS 
states in Table 1, with the CO in the ground vibrational 
state. 

Figure [3] shows the time-evolution of the population, 
P(t), for an initial wave function constructed from the 



nine zeroth order Q space states noted above, and opti- 
mized for maximal or minimal energy flow at T = 100 
fs. The optimal coefficients were found using the method 
described in Sec. [TTJ the c K coefficients and their proba- 
bilities are given in Table HH Results in panel (a) corre- 
spond to an initial superposition optimized to minimize 
the population flow from the Q to the P space, while 




0.2 



100 200 300 400 500 
t (fs) 




100 200 300 400 500 
t(fs) 



FIG. 3: IVR control in OCS: (a) IVR suppression, and (b) 
IVR enhancement. The parameters defining the optimal su- 
perposition for T — 100 fs are given in Table UU 
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panel (b) shows results optimized to enhance the flow of 
population. As is clearly seen, the initial falloff in panel 
(a) is much slower than that in (b). To quantify this 
decay, the initial P{t) falloff was fit to an exponentially 
decreasing function, 

Pit) = Poo + (I - Poo)e- t/ts , (36) 

where tg is the decay time, and is the average around 
which Pit) fluctuates for the first 1.0 ps. Note that the 
tg values can only be regarded as approximate since the 
falloff is, in general, not exponential, and tg depends on 
the time scale over which the exponential is fit. (Here 
the fit is over 400 fs). In case (a), the decay time is 
tg ~ 57.35 fs, while in case (b) it is tg ~ 8.60 fs, about 
seven times smaller. Furthermore, we note that in panel 
(a), only about 24% of the population has been trans- 
ferred from Q to P during the first 50 fs, while, in con- 
trast, approximately 82% of the population has being 
transferred to the P in panel (b) during the same time. 
Moreover, the population that asymptotically remains lo- 
calized along the CS bond is also larger in the case of 
IVR suppression (P^ ~ 0.4) than in that of enhance- 
ment (Pqo ~ 0.3). 



The controlled results should be compared to the nat- 
ural IVR behavior of the individual levels participating 
in the superposition. To this end, the P(t) for each of 
the participating levels is shown in Fig. [4j Although the 
energy difference between these states shown is relatively 
small, the populations, P K , evolve with a range of initial 
falloff values, as can be seen in the corresponding val- 
ues of tg, given in Table HT1 Note also, from this table, 
that the control seen in Fig. [3] is not due to the identifi- 
cation of a particular |k) that independently maximizes 
or minimizes the decay. Indeed, by inspecting the value 
of the c K coefficients, we find, in the case of IVR sup- 
pression, participation of most of the nine levels, with ks 
60% of the total initial population concentrated in the 
two states with n = 4 and k = 8. Neither of these two 
states independently have the longest decay times, but 
their interference is crucial to control. Similar observa- 
tions result from considering the data for optimized IVR 
enhancement, despite the fact that k = 2 has a relatively 
small tg. In this case the optimized superposition also 
gives a significantly smaller P(T) than does the individ- 
ual k = 2 state. 

A qualitative measure P(t) of the contribution from 
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FIG. 5: Contribution of P(t) (dashed line) and W(t) (dotted 
line) to: (a) IVR suppression, and (b) IVR enhancement. 
The solid line represents the corresponding P(t) function from 
Figs. El 

the interference of overlapping resonances, and W(t) 
from the direct contribution, was provided in Eq. ([58]). 
Results for P{t) and W(t) for the maximization and min- 
imization cases above are provided in Fig. [5] where the 
contribution from overlapping resonances (dashed line), 
become dominant after the first 10 fs, thus demonstrating 
the important role played by these resonances in the IVR 
control scenario. This is seen to be the case for both the 
maximization, as well as minimization, of the flow from 
the CS bond. 

A pictorial, and enlightening, view of the results is pro- 
vided in Figs. [5] and [3 where the wave packets associated 
with IVR suppression and enhancement are shown. As 
can be seen in Fig. El for the case of IVR suppression, 
the wave packet remains highly localized along the Rcs 
mode, with minimum spreading along the Rco mode. In 
particular, it undergoes a slight oscillation along the Res 
mode, concentrating most of the probability around the 
region where the CS dissociation takes place, in a clear 
correspondence to what happens with a classical counter- 
part. For the case of IVR enhancement, the effect is the 
opposite. As can be seen in Fig. [71 the spreading of the 
wave packet along the Rco mode coordinate is relatively 
fast. 

The method described above is, of course, applicable 
at any time during the dynamics. For example, we tried, 
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FIG. 6: Wave packet evolution corresponding to IVR suppres- 
sion. Dashed lines represent equipotential energy contours, 
with the innermost corresponding to the wave packet energy, 
E+ = 0.09849 a.u. 



and successfully attained, control for times at long as 
1.5 ps (corresponding to over 50 CS vibrational periods), 
resulting in about a 55% of the population localized in 
the CS bond for IVR suppression, and about 22% for 
IVR enhancement. 



V. COMMENTS AND SUMMARY 

In this paper, a method for controlling intramolecu- 
lar vibrational redistribution has been developed and has 
been applied to OCS, where extensive control over IVR 
is attained. Of particular interest is that the control is 
achieved even though the associated classical dynamics 
is chaotic. The method, wherein the coefficients of an 
initial superposition of zeroth order states are optimized, 
is shown to rely upon the presence of overlapping reso- 
nances, a feature which is expected to be ubiquitous in 
realistic molecular systems. 

We have assumed throughout this paper that the ini- 
tial state that optimizes the intramolecular vibrational 
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space, respectively, and Nt = N K + N/3. The probability 
of being in the Q space, P(t), is given by Eq. (jHJ). In order 
to find P(t), two sets of values are needed: the set of 
eigenvalues {E 7 }, and the overlap integrals a Kj7 between 
the zeroth-order states in Q and the exact eigenstates I7). 
The partitioning algorithm described below is ingenious 
in the sense that it allows one to concentrate specifically 
on obtaining these two sets of values. The method is well 
suited to small systems. 

Beginning with Eq. ifTTj) . and using Eqs. (|2"Tj) , the al- 
gorithm is as follows: 

1. Choose a starting energy i? 7 =0 , with i correspond- 
ing to the ith iteration. In particular, one may 
choose an energy close to the zeroth-order energies. 

2. Take E* from the last iteration, and compute 

3. Diagonalize rl(E^), and select one of its eigenvalues 
to be the next trial energy, E^ +1 . 



4. If 



El\ ¥ 0, go back to step 2. 



5. If \E^ +1 — E* I = 0, El^ 1 becomes the eigenvalues 
E 7 , and its corresponding eigenvector, |D 7 ), is pro- 
portional to QI7). 

6. Repeat steps 1-5 until all Nt unique eigenvalues 
E~ are obtained. 



FIG. 7: Wave packet evolution corresponding to IVR en- 
hancement. Dashed lines represent equipotential energy con- 
tours, with the innermost corresponding to the wave packet 
energy, E- = 0.09743 a.u. 



redistribution can be prepared, for a real molecule, using 
modern pulse shaping techniques. Computations display- 
ing the resultant field were not, however, carried out on 
this OCS model since they are best done using more re- 
alistic molecular potentials in higher dimensions, yield- 
ing realistic optimizing fields. Work of this kind is in 
progress. 
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APPENDIX A: NUMERICAL 
IMPLEMENTATION 



Here, we provide a route to compute the eigenvalues 
and overlap integrals via Eq. (fT7|) . We start by defining 
N K and Np to be the basis-set dimensions in the Q and P 



In the process of diagonalizing the effective Hamilto- 
nian, 7i, each eigenvector |-D 7 ) has been normalized to 
1. Therefore, the use of the algorithm leads to a loss of 
information about QI7). This makes necessary to also 
compute the constant of proportionality between QI7} 
and \Dy). This is done by requiring that (7I7) = 1 for 
the full eigenvectors. Thus, one can assert that 



Q| 7 > = C 7 |D 7 ), 



(Al) 



with C 7 being the proportionality constant. The problem 
then reduces to finding the C 7 associated with each £L. 
This is accomplished by expressing (7I7) as 

( 7 | 7 ) = ( 7 |Q|7> + (7l^l7> 

= (7|Q 2 |7> + (71^17), 

where 

(7|Q 2 |7) = |C 7 | 2 <£ 7 |i} 7 } = |C 7 | 2 , 
and, using Eq. (fT5|) . 

( 7 |P 2 | 7 ) = (jIQHPIE^-PHP}- 1 

x [E-y - PHPy 1 PHQYy). 



(A2) 
(A3) 

(A4) 



The application of the spectral resolution of an operator, 
Eq. (fT9"|). to Eq. (|A4|) leads to 



(7|P 2 |7> 



y, (j\QH\0mHQ\j) 




(e i -Ep^j 



(A5) 
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whereby, by making use of Eq. (|A1[) . one obtains 

(7|P 2 |7) = |C- y ] 2 E ( ^™ ( ! l T' y) - (A6) 



^7 — 



Now {^\P 2 \^i) is easily computed by realizing that 
(D 7 \H\f3) = Y. D l~M\H\P) 

K 

= ^D* 7 F( K |/3). 



(A7) 
(A8) 



from which one obtains the proportionality factor | C 7 1 . 

According to the procedure previously described, we 
can determine Q\j), given |-D 7 ), with the exception of a 
constant phase factor. Note that, in general, each pro- 
portionality factor, C 7 , can be written as |C 7 |e e T, where 
6 1 is a random phase. However, this is not a problem 
since the results are independent of any constant phase 
factor; as seen from Eq. ([7]), all overlap integrals appear 
in pairs, a K / 7 a* j7 , which can be expressed as 



The substitution of Eqs. jX3]) and |M| into E Q- CM 
yields 



<7i7> = i = ic 7 r 



/ 1+E E f Wl,, (A9) 



W) 



(A10) 



1 S. A. Rice and M. Zhao, in Optical Control of Molecular 
Dynamics (Wiley, New York, 2000). 

2 M. Shapiro and P. Brumer, in Principles of the Quantum 
Control of Molecular Processes (Wiley, New York, 2003). 

3 P. S. Christopher, M. Shapiro, and P. Brumer, J. Chem. 
Phys. 123, 064313 (2005). 

4 P. S. Christopher, M. Shapiro, and P. Brumer, (to be pub- 
lished) extends the treatment in References 3 and 5 to 
twenty-four mode Pyrazine. 

5 P. S. Christopher, M. Shapiro and P. Brumer, J. Chem. 
Phys. 124, 184107 (2006). 

6 H. Feshbach, Ann. Phys. (NY.) 19, 287 (1962); 43, 410 
(1967). 

7 R.D. Levine, Quantum Mechanics of Molecular Rate Pro- 
cesses (Clarendon Press, Oxford, 1969). 

8 D. Gerbasi, Ph.D. Dissertation, University of Toronto 
(2004). 

9 E. Frishman and M. Shapiro, Phys. Rev. Lett. 87, 253001 
(2001). 

10 D. Carter and P. Brumer, J. Chem. Phys. 77, 4208 (1982). 

11 M.J. Davis, Chem. Phys. Lett. 110, 491 (1984). 

12 M.J. Davis, J. Chem. Phys. 83, 1016 (1985). 

13 R. Siebert, R. Schinke, and M. Bittererova, PCCP Com- 



munications 3, 1795 (2001). 

14 The surface of section has been computed in the standard 
way, i.e., by following each trajectory, and noting Ri and 
Pi each time that R2 crosses the surface R2 = R° with 
P 2 > 0. 

15 B. Eckhardt, Phys. Rep. 163, 205 (1988). 

16 B.V. Chirikov, J. Nucl. Energy C 1, 253 (1960); Phys. Rep. 
52C, 265 (1979). 

17 R.C. Brown and R.E. Wyatt, Phys. Rev. Lett. 57, 1 (1986); 
J. Phys. Chem. 90, 3590 (1986). 

18 L.L. Gibson, G.C. Schatz, M.A. Ratner, and M.J. Davis, 
J. Chem. Phys. 86, 3263 (1987). 

19 A.J. Lichtenberg and M.A. Lieberman, Regular and 
Stochastic Motion (Springer- Verlag, Berlin, 1983). 

20 Perturbed trajectories were designed in the traditional 
manner, modifying slightly Ri (Ri —* R\ = R\ + 5) and 
adjusting P2 to restore the original energy. 

21 J.V. Lill, G.A. Parker, and J.C. Light, Chem. Phys. Lett. 
89, 483 (1982); J.C. Light, LP. Hamilton, and J.V. Lill, J. 
Chem. Phys. 82, 1400 (1985); S.E. Choi and J.C. Light, J. 
Chem. Phys. 92, 2129 (1990). 



